#!/usr/bin/env python
"""
Reads an HMM file and converts each natural log probability into a position
specific scoring matrix.

"""

from ssummo.cmd_options import Options

import math
import multiprocessing
import re
import os
import subprocess
import sys
import time

from Bio import SeqIO, AlignIO


class ColumnError( Exception ):
    def __init__(self):
        Exception.__init__(self, "\nAlignment column order wrong in HMM!!")

class LocalOptions(Options):

    options = {
            '-hmm'  : [],
            '-ncpus': multiprocessing.cpu_count() - 1,
            '-out'  : sys.stdout,
            '-in'   : None,
            '-format' : 'fasta',
            '-desc' : [],
            '-accs' : [],
            '-hmmalign' : '/biol/bin/hmmalign',
            '--consensus' : False
        }

    help_text = {
            '-hmm'   : 'HMM input file.',
            '-ncpus' : 'Number of processes to intitiate which will do the math calculations',
            '-out'   : 'Output file. Default: standard out. N.B. Appends to the file.',
            '-in'    : 'Input sequence file. Required.',
            '-desc'  : 'Filter input sequences by description(s), space delimited (Enclose with \' or " ).',
            '-accs'  : 'Filter input sequences by accession(s), space delimited.',
            '--consensus' : "Also compare each HMM consensus sequence to the model whence it came. [False]",
            '-hmmalign' : 'hmmalign program. Case sensitivity matters with the alignments!'
        }

    def print_help(self):
        print('Valid options:-')
        for key, val in self.help_text.iteritems():
            print('{0} : {1}'.format(key, val))

    def parse_args(self, args):
        options = self.options
        #global options
        if len(args) == 0:
            self.print_help()
            exit()
        prevargDash = False
        for i in xrange(len(args)):
            if args[i].startswith('--'):
                if args[i] in options.keys():
                    options[ args[i] ] = True
                else:
                    print('Invalid option: {0}'.format(args[i]))
                    self.print_help()
                    exit()
            elif args[i].startswith('-'):
                if args[i] in options.keys():
                    prevargDash = True
                    command = args[i]
                else:
                    print('Invalid option: {0}'.format(args[i]))
                    self.print_help()
                    exit()
            elif prevargDash:
                if command in ['-in', '-out']:
                    options.update( { command : args[i] } )
                    prevargDash = False
                elif command == '-accs':
                    options[command].append(args[i])
                elif command == '-desc':
                    options[command].append( '('+args[i] + ')')
                elif command == '-format':
                    options.update( {command : args[i] } )
                    prevargDash = False
                elif command == '-hmm':
                    options[command].append( args[i] )
            else:
                print("Don't know what to do with {0}".format(args[i]))
                exit()
        if options['-in'] is None:
            sys.stderr.write('"-in" is a required argument!\n')
            exit()
        return options

class AlignSeq( multiprocessing.Process ):

    def __init__(self, inPipe, outPipe, options):
        multiprocessing.Process.__init__(self)
        self.in_pipe = inPipe
        self.out_pipe = outPipe
        self.options = options.copy()
        self.options.update({'-acc':[], '-desc':[],
                            '-informat' : options['-format'] })
        self.options.update({'-format' : 'stockholm' })
        #self.sem = semaphore
        self.hmmalign = options['-hmmalign']
        self.hmmsearch = 'hmmsearch'

    def run(self):
        inp = self.in_pipe.recv()
        while inp != 'STOP':
            seq_record, HMM = inp
#            align_process = subprocess.Popen( [ self.hmmalign,'--informat',self.options['-informat'],'--outformat','A2M',HMM,'/dev/stdin'],shell=False,stdin=subprocess.PIPE,stdout=subprocess.PIPE )
            cmd = [self.hmmsearch, '-A', '/dev/stdout', '-o', os.devnull,
                   HMM, '/dev/stdin']
            align_process = subprocess.Popen(cmd, shell=False,
                                             stdin=subprocess.PIPE,
                                             stdout=subprocess.PIPE )
            align_process.stdin.write('>{0}\n'.format(seq_record.id))
            align_process.stdin.write(seq_record.seq.tostring().replace('-', ''))
            options.update( { '-in' : align_process.stdout } )
            align_process.stdin.close()
            for seq in AlignIO.read(align_process.stdout, 'stockholm'):
                self.out_pipe.send_bytes(seq.seq.tostring())
            align_process.wait()
            align_process.stdout.close()
            inp = self.in_pipe.recv()
        self.in_pipe.close()
        self.out_pipe.close()

class CalcProb(multiprocessing.Process):

    def __init__(self, inQ, sumQ, lock):
        multiprocessing.Process.__init__(self)
        self.inQ = inQ
        self.sumQ = sumQ
        self.lock = lock

    def run( self ):
        in_LO = self.inQ.get() ## Log Odds ratio parsed from HMM.
        while in_LO != 'STOP':
            while in_LO != 'endSEQ' and in_LO != 'STOP':
                try:
                    res, ResScores = in_LO.items()[0]
                    NSeqs = len( ResScores)
                    if NSeqs == 1:
                        Prob = math.exp( - float(ResScores[0]) )
                    elif NSeqs > 1:
                        Prob = math.fsum([math.exp( - float(p) ) for p in ResScores ] ) / float( NSeqs ) ## Average probabilities
                    LogOdds = - math.log( 1 - Prob ) ## Log Odds probability of residue not being observed.
                    self.sumQ.put( {res : LogOdds } )
                    in_LO = self.inQ.get()
                except Exception:
                    sys.stderr.write( str(in_LO ))
                    sys.stderr.flush()
                    self.sumQ.put( 'endSEQ' )
                    self.sumQ.put('STOP' )
                    raise
            if in_LO == 'endSEQ':
                time.sleep(0.01)
                self.sumQ.put( in_LO ) ## Send 'endSEQ'
                in_LO = self.inQ.get()

class SumProbs(multiprocessing.Process):
    def __init__(self, inQ, outQ):
        multiprocessing.Process.__init__(self)
        self.inQ = inQ
        self.outQ = outQ

    def run(self):
        in_value = self.inQ.get()
        while in_value != 'STOP':
            tally = float(0)
            probs = {}
            while in_value != 'endSEQ' and in_value != 'STOP':
                res, LogOdds = in_value.items()[0]
                tally += LogOdds   ## Sum the Log Odds probabilities.
                if res in probs.keys():
                    probs[res].append(LogOdds)
                else:
                    probs.update( {res : [LogOdds]} )
                in_value = self.inQ.get()
            if in_value == 'STOP':
                break
            else:
                self.outQ.put( [tally,probs] )
                in_value = self.inQ.get()
        self.inQ.close()
        self.outQ.close()

def parseHMM(options, prob_calc_procs):
    if options['-hmm'] == None:
        print('No HMM specified')
        for proc in prob_calc_procs:
            proc.inQ.put( 'endSEQ' )
            proc.inQ.put( 'STOP' )
            proc.join()
        exit()
    else:
        inHandle = file(options['-hmm'],'r')
    split_cols = re.compile( r'\s+' )
    line = inHandle.readline()
    HMMLineReg = re.compile(r'HMM\s+') ## Finds line that starts 'HMM' followed by white space.
    while not HMMLineReg.search(line): ## Go to that line. (Contains column headers)
        line = inHandle.readline()
    i = 0
    matchCols = {}
    for residue in split_cols.split( line.rstrip() )[1:]:  ## Column headers. For residue in header line:
        matchCols.update( { residue : i } )           #### Add { residue : column# } to matchCols.
        i += 1
    inHandle.readline()        ## Line underneath HMM line. Ignore.
    globalProbs = split_cols.split( inHandle.readline().strip() )[1:]   ## The COMPO line = Overall probabilities per residue.
    insertProbs = [split_cols.split( inHandle.readline().strip() )]        ## Insert probabilities at position 0 (before HMM).
                    #7 entries. Transition from begin / insert state 0
    begins = split_cols.split( inHandle.readline().strip() )
                    # B-->M1 ; B-->I0 ; B-->D1 ; I0-->M1 ; I0-->I0; P=1 ; P=0.
                    # ..Must use one.
    line = inHandle.readline()    # First proper line.
    con_col = 1  # Counter for consensus columns. HMMer asserts column order with first entry in each group of 3 lines.
    matchProbs = [ ]
    transitionProbs = [ ]
    while line != '//':
        entries = split_cols.split( line.strip() )
        try:
            assert int(entries[0]) == con_col
        except AssertionError:
            raise ColumnError(entries[0])
        matchProbs.append( entries[1:] ) ## Making the array. Hiccup: indexes are 1 No. below column No. 1 per residue & Col# ; RF Annotation ; Consensus structure.
        insertProbs.append( split_cols.split( inHandle.readline().strip() ) ) # 1 per residue.
        transitionProbs.append( split_cols.split( inHandle.readline().strip() ) ) # Mn-->Mn+1 ; In ; Dn+1 ; In-->Mn+1 ; In ; Dn-->Mn+1 ; Dn+1
        try:
            assert len( transitionProbs[-1] ) == 7
        except Exception:
            sys.stderr.write( 'Transition line (below) does not have 7 states...' )
            sys.stderr.write( transitionProbs[-1] )
            raise
        line = inHandle.readline().strip()
        con_col += 1
    matrices = {
            'col_index' : matchCols,
            'match' : matchProbs,
            'insert': insertProbs,
            'globals': globalProbs,
            'begin'    : begins,
            'transitions' : transitionProbs
            }
    inHandle.close()
    return matrices

def score_alignment( alignment_residues, matrices, processes ):
    IUPACcodes = {     ##'T' : ['U'], ## Will consistently only turn U --> T.
            'T' : ['T'],
            'A' : ['A'],
            'G' : ['G'],
            'C' : ['C'],
            'U' : ['T'],
            'M' : ['A', 'C'],
            'R' : ['A', 'G'],
            'W' : ['A', 'T'],
            'S' : ['C', 'G'],
            'Y' : ['C', 'T'],
            'K' : ['G', 'T'],
            'V' : ['A', 'C', 'G'],
            'H' : ['A', 'C', 'T'],
            'D' : ['A', 'G', 'T'],
            'B' : ['C', 'G', 'T'],
            'X' : ['G', 'A', 'T', 'C'],
            'N' : ['G', 'A', 'T', 'C']
            }
    col_ind = 0
    proc_end = len(processes)
    proc_ind = 0
    residues = matrices['col_index']
    if 'U' in residues:
        residues.update( { 'T' : residues['U'] } )
    state = 'begin'
    col_ind = 0
    res_ind = 0
    seq = alignment_residues
    ## B --> M1, B --> I0, B --> D1; I0 --> M1, I0 --> I0; then a 0.0 and a .*.
    while state != 'match' and state != 'delete':  ## Look through until leave begin state.
        probs = []
        res = seq[res_ind]
        if res in ['-', '~']:        ## Delete state.
            processes[proc_ind].inQ.put( { col_ind :[matrices['begin'][2]] } )
            proc_ind += 1
            state = 'delete'
            col_ind += 1
        elif res in ['.']:    ## Insert state
            if state == 'insert':
                processes[proc_ind].inQ.put( { col_ind : matrices['begin'][4] } )
            elif state == 'begin':
                processes[proc_ind].inQ.put( { col_ind : matrices['begin'][1] } )
            state = 'insert'
            proc_ind += 1
        elif res.isupper():    ## Match state.
            assert col_ind == 0
            if state == 'begin':
                processes[proc_ind].inQ.put( {col_ind :[matrices['begin'][0]] } )  ## B --> M1
            elif state == 'insert':
                processes[proc_ind].inQ.put( {col_ind :[matrices['begin'][3]]} )  ## I0 --> M1
            else:
                raise ValueError("The state {0} is not valid for before the first match state".format(state))
            for trans in IUPACcodes[res.upper() ]:
                probs.append( matrices['match'][col_ind][ residues[ trans] ] )
            processes[proc_ind].inQ.put( { col_ind : probs } )
            state = 'match'
            col_ind += 1
            proc_ind += 1
        else:            ## Lower case == insert state.
            if state == 'insert':
                processes[proc_ind].inQ.put( { col_ind : [ matrices['begin'][4] ] }) ## I0 --> I0
            elif state == 'begin':
                processes[proc_ind].inQ.put( { col_ind : [ matrices['begin'][1] ] })     ## B-->I0
            else:
                sys.stderr.write( "Can't go from state {0} to insert".format(state) )
            state = 'insert'
            for trans in IUPACcodes[res.upper() ]:
                probs.append( matrices['match'][col_ind][ residues[trans] ] )
            processes[proc_ind].inQ.put( {col_ind :probs} )
        if proc_ind >= proc_end :
            proc_ind = 0
        res_ind += 1

    for res in seq[col_ind:]:
        try:
            probs = []
            if res in [ '-', '~'  ]: # Allowed gap residues (and probably more...).
                if state == 'match':
                    processes[proc_ind].inQ.put( {col_ind : [matrices['transitions'][col_ind][2]] } ) ## Mk --> Dk+1
                elif state == 'delete':
                    delVal = matrices['transitions'][col_ind][6]
                    if delVal != '*':
                        processes[proc_ind].inQ.put( {col_ind : [matrices['transitions'][col_ind][6]] } ) ## Dk --> Dk+1
                else:
                    sys.stderr.write( "Can't go from insert state to delete state!" )
                    sys.stderr.flush()
                state = 'delete'
                col_ind += 1
            elif res == '.':      # Insert state
                if state == 'insert':
                    in_prob = matrices['transitions'][col_ind][4]
                elif state == 'match':
                    in_prob = matrices['transitions'][col_ind][1]
                else:
                    sys.stderr.write( "Can't go from state {0} to 'insert'".format(state) )
                processes[proc_ind].inQ.put( {col_ind : [in_prob] } )
                for trans in IUPACcodes[res.upper() ]:
                    probs.append( matrices['insert'][col_ind][ residues[trans] ] )
                state = 'insert'
            elif res.isupper():
                if state == 'match':
                    processes[proc_ind].inQ.put( {col_ind : [matrices['transitions'][col_ind][0]]} ) ## Mk --> Mk+1
                elif state == 'delete':
                    processes[proc_ind].inQ.put( {col_ind : [matrices['transitions'][col_ind][5]]} ) ## Dk --> Mk+1
                elif state == 'insert':
                    processes[proc_ind].inQ.put( {col_ind : [matrices['transitions'][col_ind][3]]} ) ## Ik --> Mk+1
                else:
                    sys.stderr.write( "Unknown state: {0}.\n".format(state) )
                for trans in IUPACcodes[ res ]:
                    probs.append( matrices['match'][col_ind][ residues[trans ] ])
                state = 'match'
                col_ind += 1
            else:
                if state == 'insert':
                    in_prob = matrices['transitions'][col_ind][4]  ## Ik --> Ik
                elif state == 'match':
                    in_prob =  matrices['transitions'][col_ind][1]  ## Mk --> Ik
                else:
                    sys.stderr.write( "Can't go from state {0} to 'insert'".format(state) )
                processes[proc_ind].inQ.put( {col_ind : [in_prob] } )
                try:
                    for trans in IUPACcodes[res.upper() ]:
                        probs.append( matrices['insert'][col_ind][ residues[trans] ] )
                except Exception:
                    print(col_ind)
                    print(len( matrices['insert'] ))
                    print(residues)
                    print(matrices['insert'][col_ind])
                    raise
                state = 'insert'
            proc_ind += 1
            if proc_ind >= proc_end:
                proc_ind = 0
            processes[proc_ind].inQ.put( {col_ind :probs } ) ## -ve log odds probability.
            res_ind += 1
        except Exception:
            print('\ncol ind', col_ind)
            print('len insert: ', len(matrices['insert']))
            print('len matches: ', len(matrices['match'] ))
            print('len transi: ', len(matrices['transitions']))
            print('seqlen: ', len(seq))
            print(seq)
            print('residues', residues)
            print('matrices[insert][col_ind]:,', matrices['insert'][col_ind])
            #print 'matrices[match][col_ind]:,', matrices['match'][col_ind]
            print('proc_ind:', proc_ind)
            print('proc_end:', proc_end)
            raise
    processes[proc_ind].inQ.put( 'endSEQ' )
    return True


def get_consensus( HMM ):
    emit_process = subprocess.Popen( [ 'hmmemit', '-c', HMM ], shell = False, stdout = subprocess.PIPE )
    consensus_seq = SeqIO.read( emit_process.stdout, 'fasta' )
    emit_process.stdout.close()
    return consensus_seq

def parse_seqs(options):
    if type(options['-in']) != file:
        options.update( { '-in' : file(options['-in'], 'r') } )
    if len(options['-desc']) > 0:
        desc_reg = re.compile( '|'.join([arg for arg in options['-desc']] ), re.I )
    noFilter = False
    if len(options['-accs']) > 0:
        acc_reg = re.compile( '|'.join([re.escape(arg) for arg in  options['-accs']] ), re.I )
    elif len(options['-accs']) == 0 and len(options['-desc']) == 0:
        noFilter = True
    Ndesc = len(options['-desc'])
    Naccs = len(options['-accs'])
    count = 0
    for seq in SeqIO.parse( options['-in'], options['-format']  ):
        count += 1
        if Naccs > 0:
            accFind = acc_reg.search( seq.id )
            if accFind:
                sys.stdout.flush()
                yield seq
                continue
        if Ndesc > 0:
            if desc_reg.search( seq.description ):
                yield seq
        elif noFilter:
            yield seq

def print_seqs_vs_HMM(accession_names, matrices, out=sys.stdout):
    """This prints multiple seqeunces against each sequential HMM. So the output
    will contain columns containing the cumulative sequence scores of each sequence
    against just one model. If there's multiple HMMs, then the next HMM will be
    printed underneath.
    """
    HMMs = sorted( matrices.keys() )
    accessions = sorted(accession_names.keys())
    out.write('\n')
    out.write('# Each section is split by headers where each line represents a different descriptive header.\n')
    out.write('# The first line gives the name of the HMM.\n')
    out.write('# The second line always the accession numbers for the target sequences.\n')
    out.write('# The next line will give the first 30 characters in the descriptive sequence header.\n')
    out.write('# After that comes all the marvelous data....\n')
    out.write('\n# These HMMs will be searched, in this order: ')
    out.write( str(', '* len(accessions)).join(  HMMs ) +'\n\n')
    if 'cons' in accessions:
        del(accessions[ accessions.index('cons') ] )
        accessions.append('cons')
    for n in xrange( len(HMMs) ):
        ## Print column headers
        out.write( '# ***|' + str('\t'*len(accessions)) + HMMs[n] + str('\t'*len(accessions))+'|***\n')
        out.write( '#\t'+ '\t'.join( [accession_names[acc][:30] for acc in accessions ] ) +'\n')  # Sequence description line
        out.write('#Col No.\t'+ ( '\t'.join(accessions))+'\n')                    ## Sequence accession
        tallies = [ float(0) for i in xrange( len(accessions ) )]
        matrix = matrices[HMMs[n]]
        print([len(matrix[acc][1].keys()) for acc in accessions])
        lengths = [ len(matrix[acc][1].keys()) for acc in accessions ]
        nResidues = max( lengths )
        for i in xrange( nResidues):
            out.write('\n{0}\t'.format(i))
            #colInds = {}
            #inserts = []
            for seqInd in xrange(len(accessions)):
                if lengths[seqInd] <= i:
                    matrix[ accessions[seqInd] ][1].update( { i: [] } )
                for score in matrix[accessions[seqInd]][1][i]:
                    tallies[seqInd]  += score
            out.write('\t'.join([str(t) for t in tallies]))
        out.write('\n\n')

def print_HMMs_vs_seq(accession_names, matrices, resultseq, out=sys.stdout):
    """This will write (to out) the cumulative sequence scores a single sequence against
    multiple HMMs. The columns are subdivided into groups, each group corresponding to an
    HMM. Within each group of columns, every sequence gets a column of it's own.
    Note that consensus sequence scores show the cumulative scores for consensus sequences
    that are unique for each HMM."""
    HMMs = sorted( matrices.keys())
    accessions = sorted(accession_names.keys())
    if 'cons' in accessions:
        del(accessions[ accessions.index('cons') ] )
        accessions.append('cons')
    out.write('# Each section is split by headers where each line represents a different descriptive header.\n')
    out.write('# The first line gives the name of the HMM.\n')
    out.write('# The second line always the accession numbers for the target sequences.\n')
    out.write('# The next line will give the first 30 characters in the descriptive sequence header.\n')
    out.write('# After that comes all the marvelous data....\n')
    out.write('\n# These HMMs will be searched, in this order: ')
    out.write( ', '.join( HMMs ) +'\n')  # Sequence description line
    out.write( '\t'.join(accessions) +'\n\n')
    if 'cons' in accessions:
        del(accessions[ accessions.index('cons') ] )
        accessions.append('cons')
    nSeqs = len(accessions)
    HMMLens = []
    maxHMMLen = 0
    for HMM in HMMs:
        HMMLens.append(  [len(matrices[HMM][accession][1].keys()) for accession in accessions]  )
        maxHMMLen = max( [maxHMMLen] + HMMLens[-1] )
    #HMMLens = [ len(matrices[HMM][accessions[0]][1].keys()) for HMM in HMMs]
    nHMMs = len(HMMs)
    ResInd = 0
    tallies = [ [float(0) for i in xrange(nSeqs)] for HMM in xrange(len(HMMs))]
    tallyInd = 0
    #insertReg = re.compile( r'[a-z\.]{1}' )
    try:
        while ResInd < maxHMMLen:
            line = []
            HMMInd = 0
            while HMMInd < nHMMs:
                SeqInd = 0
                while SeqInd < nSeqs: # Each sequence gets a column in each row.
                    if ResInd < HMMLens[HMMInd][SeqInd]:
                        scores = matrices[ HMMs[HMMInd] ][ accessions[SeqInd] ][1][ ResInd ]
                        tallies[HMMInd][SeqInd] += sum( scores )
                        line.append( tallies[HMMInd][SeqInd] )
                    else:
                        line.append( tallies[HMMInd][SeqInd] )
                    SeqInd += 1
                    tallyInd = SeqInd
                HMMInd += 1
            out.write('{0}\t'.format(ResInd)+ '\t'.join( str(score) for score in line ) +'\n')
#            if insertReg.search( resultseq[ResInd] ):
#                continue
#            else:
            ResInd += 1
    except Exception:
        print(HMMLens, HMMs, HMMInd, SeqInd, nSeqs,
              accession_names.keys(), matrices.keys())
        print(accessions[SeqInd])
        print("ResInd: ", ResInd)
        print("MaxHMM Len: ", HMMLens[HMMInd])
        raise
        count = 0
        err = sys.stderr.write
        for i in sorted(matrices[HMMs[HMMInd]][accessions[SeqInd]][1]):
            try:
                assert i == count
            except AssertionError:
                err('i: {0}, count: {1}, len: {2}\n'.format( i, count,
                    len(matrices[HMMs[HMMInd]][accessions[SeqInd]][1])))
                raise
            count += 1
        print("N Consensus Cols: ", count)
        raise


if __name__ == '__main__':
    args = sys.argv[1:]
    options = LocalOptions().parse_args( args )
    #### PROCESS INITIATION ####
    ## Start Qs & Pipes between SumProbs & CalcProb ####
    calc_to_sumQ = multiprocessing.Queue()
#    sum_end, calc_end = multiprocessing.Pipe()
    sum_seq_queue = multiprocessing.Queue()
    ## Initiate the SumProbs process ##
    sum_thread = SumProbs( calc_to_sumQ, sum_seq_queue )
    sum_thread.start()
    lock = multiprocessing.RLock()
    ## Initiate the CalcProb processes ##
    math_processes = [ CalcProb(multiprocessing.Queue(), calc_to_sumQ, lock) \
                        for i in xrange(options['-ncpus']) ]
    ## hmmalign controller processes & semantics ##
    seq_send_pipe, seq_recv_pipe = multiprocessing.Pipe()
    ali_recv_pipe, ali_send_pipe = multiprocessing.Pipe()
#    semaphore = multiprocessing.Semaphore(options['-ncpus'])
#    align_processes = [ AlignSeq( seq_recv_pipe, ali_send_pipe, semaphore, options) for i in xrange(options['-ncpus'] ) ]
    align_process =  AlignSeq( seq_recv_pipe, ali_send_pipe, options)
    for i in xrange(options['-ncpus']):  ## Start the bad boiiis
        math_processes[i].start()
    align_process.start()
    ####  END PROCESS INITIATION #####
    sequence_scores  = {}
    acc_names = {}
    sys.stdout.write( '# '+ '\t'.join( ['Accession', 'HMM', 'Score'] ) )
    if options['--consensus']:
        acc_names.update( { 'cons' : 'HMM Consensus sequence' } )
        sys.stdout.write( '\t{0}'.format( 'Cons. Score' ) )
    env = options.copy()
    for seq in parse_seqs( options ):  ## Iterate through sequences
        acc_names.update( { seq.id : seq.description[len(seq.id):].strip() } )
        for HMM in options['-hmm']:    ## Iterate through HMMs.
            env.update( { '-hmm' : HMM } )
            matrices = parseHMM( env, math_processes )  ## Convert the HMM to log odds probability matrices.
            seq_send_pipe.send([seq, HMM])           ## Send sequence to be aligned.
            ali = ali_recv_pipe.recv_bytes()       ## Receive aligned seqeunce.
            submitted = score_alignment( ali, matrices, math_processes )  ## Score sequence against HMM.
            tally, probs = sum_thread.outQ.get()       ## Get the results.
            HMMFname = HMM[HMM.rfind(os.path.sep)+1:]  ## Don't need the full path name; just the file name when writing out.
            if HMMFname in sequence_scores.keys():
                sequence_scores[HMMFname].update( { seq.id : [tally, probs] } )
            else:
                sequence_scores.update( { HMMFname : { seq.id : [tally, probs] } } )  ## Update dictionary of results.
            sys.stdout.write('\n# {0}\t{1}\t{2}'.format( seq.id, HMMFname, tally ) )
            if options['--consensus']:    ## Do same for the consensus sequence (if '--consensus' )
                if 'cons' not in sequence_scores[HMMFname].keys():
                    cons_seq = get_consensus( HMM )
                    submitted = score_alignment( cons_seq, matrices, math_processes )
                    cons_tally, cons_probs = sum_thread.outQ.get()
                    sequence_scores[HMMFname].update( {'cons' : [cons_tally, cons_probs]} )
                sys.stdout.write('\t{0}'.format( sequence_scores[HMMFname]['cons'][0] ) )
            sys.stdout.write('\t{0}'.format( len( seq )  ))
            sys.stdout.flush()
    ##### KILL PROCESSES #####
    for i in xrange( options['-ncpus'] ):
        math_processes[i].inQ.put( 'STOP' )
    seq_send_pipe.send( 'STOP' )
    align_process.join()
    for i in xrange( options['-ncpus'] ):
        math_processes[i].join()
    calc_to_sumQ.put( 'STOP' )
    sum_thread.join()
    ##### PRINT MATRIX ####
    if options['-out'] == sys.stdout:
        save = raw_input("\nProcessed all sequences. Do you want to print the distance matrices? [y, n or filename]" )
        if re.search(r'^ye?s?$', save, re.I):
            print_HMMs_vs_seq( acc_names, sequence_scores, ali, out=sys.stdout )
        elif re.search(r'^no?$', save, re.I) or save.strip() == '':
            exit()
        else:
            sys.stdout.write( 'Saving tab-delimited data to {0}.\n'.format( save ) )
            outFile = file(save, 'w')
            print_HMMs_vs_seq(acc_names, sequence_scores, ali, out=outFile)
            outFile.close()
    else:
        sys.stdout.write('\nSaving tab-delimited data to {0}.\n'
                         .format(options['-out']))
        if os.path.exists( options['-out'] ):
            sys.stdout.write('Overwriting {0}.\n'.format(options['-out']))
        outFile = file(options['-out'], 'w')
        print_HMMs_vs_seq( acc_names, sequence_scores, ali, out=outFile )
    options['-in'].close() ## At end in case this refers to standard input, in which case raw_input won't work.
